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Abstract 

In this article, we extend a Milstein finite difference scheme intro- 
duced in Giles Reisinger(2011)] for a certain linear stochastic par- 



tial differential equation (SPDE), to semi- and fully implicit timestep- 
ping as introduced by [Szpruch(2010)| for SDEs. We combine stan- 
dard finite difference Fourier analysis for PDEs with the linear stability 
analysis in Buckwar Sz Sickenberger(2011)| for SDEs, to analyse the 



stability and accuracy. The results show that Crank-Nicolson time- 
stepping for the principal part of the drift with a partially implicit 
but negatively weighted double ltd integral gives unconditional stabil- 
ity over all parameter values, and converges with the expected order in 
the mean-square sense. This opens up the possibility of local mesh re- 
finement in the spatial domain, and we show experimentally that this 
can be beneficial in the presence of reduced regularity at boundaries. 

Keywords: Stochastic partial differential equations, finite differences, implicit 
timestepping schemes, Fourier analysis, local mesh refinement 



1 Introduction 

The numerical analysis and computation of stochastic partial differential 
equations (SPDEs) have become a subject of active research over the re- 
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cent past. The interest has been triggered partly by applications in ar- 
eas as divers e as geophysics [Winter &: Tartakov sky(200 2~)] and m athemat- 
ical finance Heath et al. (1992), Musiela fc Zariphopoulou(2009)| , and has 
led to questions regarding the complexity theory of their approximation 
Miiller-Gronbach & Ritter(2007)[ [Miiller-Gronbach et a/.(2007)| . 

In a prominent class of SPDEs, the stochasticity enters via a random 



driver of the form 



dv = {Av + f(v))dt + g(v) dM t , 



(1) 



where A is a linear elliptic operator, e.g. the Laplace operator, M a martin- 
gale driver, often standard Brownian motion, and / and g non-linear func- 
tions, e.g. with Lipschitz regularity. This leads to special cases with additive 
or multiplicative noise terms. We will consider here a variant of equation ([1]) . 

Typical solutions are by lattice methods, e.g. Gyongy(1999) , Gyongy fc Nualart(1997)| , 
finite differences, e.g. |Roth(2002)| , or by finite ele ments, see e.g. |Walsh(2005)| , 



with extensions to higher order Taylor schemes Jentzen & Kloeden(2009) 



Jentzen & Kloeden(2010) Jentzen et a/. (2011)] , as well as multilevel schemes 
Barth et a/.(2011)] . 



In the following, let (fi, J 7 , P) a probability space, M a one- dimensional 
standard Brownian motion adapted to J 7 . We study specifically the equation 



dv , 1 d 2 v , _ dv , , , 



(2) 



where \x and < p < 1 are real-valued parameters. It is a classical re- 
sult from Krylov fc Rozovskii(198i)] that for a class of parabolic SPDEs 
including (jjj), with initial data in H 1 , there is a unique weak solution v G 
L 2 {VL x (0,T),J T ,H 1 (R)). In fact, for the special form © on R, i.e. without 
boundaries, it is easy to see that a solution is given by 

v(t, x) = u(t, x — fit — y/pM t ), 



where u is the solution to the heat equation 



du 
dt 



p) 



d 2 u 
dx 2 



with the same initial data as (F2J). We will use this semi-analytical solution 
to measure the errors of numerical approximations. 

We list two applications of this equation. Kurtz &: Xiong(1999)| show 
that 02]) describes the limit empirical measure of a large particle system, 
where each individual particle is governed by 



dX l t = fx dt + y/1 - p dW l t + yfp&Mt 



(3) 
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where W % are standard Brownian motions, which are independent mutually 
and of the Brownian motion M. The parameter p describes the correla- 
tion between each pair of X 1 , which explains the motivation for choosing 
< p < 1 in (j2J). We will see later that p = 1 is a boundary case in the sta- 
bility analysis. It is also clear how equations with non-normalised constant 
coefficients can be rescaled to fl2]). 

Equation d2J) also arises as the Zakai equation in a stochastic filtering 
problem (see, e.g. Bain fc Crisan(2009)| ), where the solution is the distribu- 
tion of a signal X, based on noisy observation of M. 

Giles &: Reisinger(2011)| introduce a Milstein finite difference approxi- 



mation of (jzj) and study the complexity of multi-level Monte Carlo simulation. 
In this article, we extend the discretisation and its analysis to an implicit 
method on the basis of the a-6 time-stepping scheme proposed and analysed 
by [Szpruch(2010) Buckwar fc Sickenberger(2011)| for SDEs, where the drift 
and the deterministic part of the double stochastic integral are taken (partly) 
implicit. By combining Fourier methods in Buckwar fc Sickenberger(2011)| 
and |Giles fc Reisinger(2011)| , we obtain a stability condition on the ratio 
k/h 2 , under which the approximations to the initial- value problem of ([2]) con- 
verge in mean-square sense in the spatial L%- and Loo-norms, of first order 
in the time-step k and second order in the spatial mesh size h. A peculiarity 
of equation (T2]) is that the stability region of the chosen scheme is larger for 
explicit treatment of the Milstein correction than for partly implicit treat- 
ment as explained above, and that an 'anti'-implicit version with negative 
weight of the implicit term gives unconditional stability. We find, both from 
an asymptotic expansion of the error and numerical experiments, that the 
numerical error is dominated by the stochastic terms of the equation and 
therefore implicit or even Crank-Nicolson-type versions of the scheme have 
little effect on the achieved accuracy. The ratio of k/h 2 , which gives, empiri- 
cally, optimal accuracy for constant mesh sizes, is close to the stability limit 
of the explicit scheme. The improved stability can, however, be useful for 
locally refined schemes. 

The rest of the article is organised as follows. In Section (2J we define 
the implicit Milstein finite difference schemes and analyse their stability and 
accuracy by Fourier techniques. Section[3]presents numerical tests which con- 
firm and illustrate these findings. Section H] gives an application to the pric- 
ing of basket credit derivatives, where the presence of an absorbing boundary 
leads to local loss of regularity, and we show how mesh grading in conjunc- 
tion with unconditionally stable implicit schemes can be used to improve 
the accuracy. Section \5\ gives conclusions and outlines directions for future 
research. 
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2 Discretisation and analysis of stability and 
convergence 

Starting point is the integrated form of the SPDE (J2j), over a time interval 

[t,t+k], 



f t+h ( dv 1 d 2 v\ , f t+k ^dv 



2.1 Milstein finite differences 

In Giles fc Reisinger(2011)| , a Milstein approximation to the stochastic in- 
tegral is used, together with standard central difference approximations on 
a spatial grid with uniform spacing h, to obtain an approximation vj to 
v(nk,jh) defined by 



where Z n ~ iV(0, 1) are independent, for n > 0. 

For a vector V n = (. . . , u ™ 1; Uq, t> v% , ■ ■ •) € K z , the system can then be 
written in operator form 

t/ t/ yk + y/pkZ n {l-p)k + pkZl 

K+i = K - DiV n + — 2 D 2 V n , (5) 

where D\ and D 2 are first and second central difference operators. 

Remark 2.1 The discretisation arises from a 'horizontal' method of lines, 
where the time integral is approximated first, and then the spatial derivatives 
are approximated by finite differences. The 'vertical' version where the Mil- 
stein scheme is applied to the system of SDEs resulting from a finite difference 
approximation of the spatial derivatives, leads to 

p:k + y/pkZ n k pk(Zl-l) 2 



V n+l = V. - - " D lV . + W2 D,V. + ^j^D-X. (6) 

The only difference is in the ltd term, where the second difference is replaced 
by an iterated first difference. We will sketch in Remark \2.2\ why the prop- 
erties of the schemes are asymptotically identical, while the scheme proposed 
earlier has implementational advantages as it leads to more compact finite 
difference stencils. 
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Giles fc Reisinger(2011)| derive the condition (l+2p 2 )k/h 2 < 1 for mean- 



square stability of this explicit scheme. For p — 0, this reduces to the well 
known stability condition for the standard heat equation. The limitation on 
the timestep is the motivation for considering the following implicit versions. 
In the spirit of K loeden &: Platen(1992)| , pp. 399, we define an implicit 



Milstein finite difference scheme by 

K+i = Vn ■ ^-D 1 V n+1 + ^D 2 V n 
2h 2h z 



2h 1 " 2h 2 

All drift terms are taken implicit, while the volatility terms are taken explicit. 
We also define a scheme 

Vn+1 = K-0 - ±D,) V M -(1-0) gl,, - ±D 2 ) V, 

^. DlVn+ eh&zllw., (8, 



2h 1 " 2h 2 

for 9 G [0, 1]. Clearly, for = one recovers the explicit scheme, for 6 = 1 
the implicit scheme. 

It is pointed out in Higham(2000b)| that the stability region of drift- 



implicit Milstein schemes is often lower than their Euler-Maruyama coun- 
terpart. |Szpruch(2010)| and Buckwar fc Sickenberger(2011)| suggest a a-9- 



scheme, which translates into the present SPDE setting as 

v n+ i = v n - 9 ( ^ a - JLd*) v n+1 -{i-o)( ^-d x - Ad, ) \ 



v 2h 1 2h 2 7 y ' \ 2h 1 2h 2 

-a§D 2 V n+1 -(l-a)^D 2 V n 

- y^D lVn + ^D 2 V n , (9) 
2h 2h 2 ' v ; 

where the deterministic part of the double Ito integral is also taken partly 
implicit with a G [0, 1]. Note that all terms that can be taken implicit, 
consistent with the Ito integral, are taken implicit. Implicitness of terms 
involving M changes the character of the integral, e.g. in the Stratonovic 
sense for a trapezium rule approximation. 

We now analyse accuracy and stability of the above schemes. The analysis 
is done on the real line (infinite grid), for analytical tractability, although in 
practical applications truncation to a finite domain and approximation on a 
finite grid will be necessary. We outline this at the start of Section |3] and 
discuss the boundary behaviour in Section HI 
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2.2 Mean-square stability analysis of Fourier modes 

We assume for simplicity fi = in the following, but the results are unaltered 
in the case \i 7^ 0, as we will discuss briefly in Remark 12.31 



As per classical finite difference analysis, e.g. Richtmyer & Morton(1967) 



Morton fc Mayers(2005)| , we study simple Fourier mode solutions of the form 

V j n = X n exp{ij<f>), |0| <vr. (10) 

We use superposition of these solutions for different to construct the leading 
order error terms in the next section, and for now focus on the stability of 
individual modes. 

Following |Higham(2000a)]|Saito fc Mitsui (1996)] , we say that the scheme 
is mean-square stable, if for all 0^0, for the amplitude X n of the corre- 
sponding Fourier mode, 

E [|X n | 2 ] ^ for n ->• 00. (11) 

Theorem 2.1 Assume p G [0, 1). The 9-a Milstein central difference scheme 
(TJP is stable in the mean-square sense ( TO]) for Fourier modes (G2P, provided 

^f(p;6,a):=^[l-2(6-pa-p 2 )] < 1. (12) 

Proof Insertion of (|10p in ([9]) leads to the equation 

X n+1 = X n + ka(9X n+1 + (l-6)X n )-VkicZ n X n (13) 
+ kpa Z 2 n X n - kpa (aX n+1 + (1 - a)X n ) , (14) 

where 

2 n 

a = -^sm -, (15) 

c = ^^sin0. (16) 
h 

Rearranging, taking moduli and expectations gives 

(l-ka(9-pa)) 2 E[ |X n+ i| 2 ] = ((l + ka{l-9+p<r)) 2 + kc 2 + 2k 2 p 2 a 2 ) E [ \X n \ 2 } . (17) 

Simple calculations show that the scheme is stable in the above mean-square 
sense if and only if 

l + Jfea(±-* + pff + ^)>-£. 
Re-inserting a and c shows that this is equivalent to 

h_ [1 _ 2(9 -pa- p 2 )] sin 2 | + pcos 2 ± < 1, V 0, 
leading to the result. □ 
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Remark 2.2 These results are similar but not identical to those obtained 
for scalar complex-valued test equations in Buckwar & Sickenberger(20ll)\ , 



because the discretisation (TJ|) differs from the standard Milstein scheme for 
systems of SDEs, as per Remark \2.1\ The stability conditions only differ by 
terms which vanish as k, h — )■ with k/h 2 fixed, and are hence asymptotically 
equivalent. This was to be expected given the schemes differ only in the stencil 
width for the discretisation of the second derivative in the ltd term. We do 
not reproduce the analysis of the other scheme here. 

Remark 2.3 Similarly, the inclusion of a first order term in the drift, p ^ ; 
leads to lower order corrections in k, k/h 2 fixed, and therefore does not change 
the result asymptotically. 

The scheme is unconditionally mean-square stable, i.e., without condi- 
tions on k and h, if / < 0. For all other cases the scheme is mean-square 
stable if k/h 2 < 1/f. This upper bound 1/f is shown as a function of p in 
Figure [U for the explicit, implicit and double implicit schemes. 

The stability condition ( Tl2|) is stricter for a > than for a = 0, i.e. scheme 
(J7J). Specifically, the doubly implicit Milstein scheme with 9 = a = 1, 
is unconditionally stable in the mean-square sense only if p < 1/(1 + VS), 
whereas the drift-implicit scheme (i.e., a = 0, 9 = 1) is unconditionally stable 
for p < l/y/2. 

This arises from the fact that the implict discretisation of the Ito term 
on the second line of containing D 2 V n+1 , has the opposite sign of the 
implicit D 2 V n+ \ term on the first line, which arises from the discretisation of 
the u xx term in the SPDE (T5]). The latter determines the parabolic nature of 
the problem. Hence, increasing a reduces this component in the implicit term 
while it increases it in the explicit term, making the scheme less contractive 
in the mean-square sense for all non-zero wave numbers, as is eventually seen 
from (1171) . Conversely, taking a < improves the stability, and for a = — 1, 
9 > 1/2, the scheme is unconditionally stable for all < p < 1. This 
somewhat surprising feature arises due to the purely imaginary eigenvalues 
of the first order operator in the Brownian driver. 

In the case 9 = 1/2, specifically, the time discretisation of the PDE part 
is of second order accurate and the error is thus dominated by the Milstein 
discretisation of the stochastic integral. 

2.3 Fourier analysis of mean-square convergence 

We can also derive the leading order error terms exploiting the availability of 
a closed-form solution. This is different from the approach in Lang(2010)| 
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Figure 1: Shown are the stability regions for the explicit scheme (9 = a — 0), 
implicit (9 = 1, o = 0) and double implicit scheme (9 = a = 1). The implicit 
scheme is unconditionally stable for p < l/y2 ~ 0.7, the double implict 
scheme for p < 1/(1 + \/3) 0.36, marked by vertical lines. In all other 
cases, the curve of the style defined in the legend gives the upper limit of the 
stability range of k/h 2 , i.e., 1/ f(p; 9, a) with / from (TT2j) . as a function of p. 

who shows a Lax-equivalence theorem for a different class of SPDEs to deduce 
convergence from stability and stochastic consistency. 

Theorem 2.2 Assume p G [0, 1), T > 0, k = T/N and £ = A > is kept 
fixed such that (TjJ^ holds. The 9-a Milstein central difference scheme |3J) 
has the error expansion, for Dirac initial data, 

Vf - v(T,jh) = kE(T,jh) + o(k) R(T,jh), (18) 

where E and R are random variables with bounded moments. 

Corollary 2.1 Under the conditions of Theorem \2.B. the 9-a Milstein scheme 
converges in the mean-square sense for the spatial L 2 - and L^-norms, and 

y/E[\\V»-v(T,.W] = 0(k). 

Proof [of Theorem 12. 2 j By insertion one checks that X (t) exp(iKx) , is a solution 
to © iff 

X(t) = X(0) exp (-|(1 - p)K 2 t - iKs/pMt) . 
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This allows us to compare the numerical solution from (|13p . 

1 + fcg(l - 0) - \/facZ ra - fc/3a((l - a) - Zp 
n+1 n 1 — kaO + /cpacr 

with a and c as in (I15p and (|16fl . to the exact solution over a single timestep, 
X{t n+ i) = X{t n ) exp (-1(1 - p)n 2 k - iK\fpkZ^ , 

where M tn+1 - M tn = V~k Z n . 

We extend here the analysis in |Giles fe Reisinger(2011)] for the explicit scheme, 
8 = a = 0, both in scope and in detail. In the following, we keep k/h 2 = A fixed, 
and set (ft = hn. 

We first consider the regime k < h~ m , where m < 1/2. Then x = kn 2 = Xh 2 K 2 
is small and one can derive by Taylor expansion 

1 + kaA - VkicZ n + kpaZ 2 = 1 - xY n (l - - i^[x^fpZ n {\ - ^) + o(x 2 ), 

= exp ^— \AkK 2 — i\fp~k,KZ n + e^ ^ , 

where A is a fixed constant, Y n = ^(A + pZ 2 ) and 

4°> = -ix^Tp [(Y n - ±)Z n - § Z*] + x 2 [(Y n - ±)(pZ 2 - \Y n ) - t*Z*\ + o(x 2 ). 
Similarly, for fixed B, 

1 - kaB = exp (±Bx - «4 1} ) , 

The remainder term o(x 2 ) = o(k 2 K 4 ) = o(kh 2 K 4 ) is understood to be a determin- 
istic constant of order o(x 2 ) multiplied by a random variable whose moments are 
all bounded independent of k and h. 

Taking A = (1 — 9) — p(l — a) and B = 6 — per, we have A + B = 1 — p and 

X n+ i = X n exp (— 1(1 — p)K 2 k — iK\J~pkZ n + e 



where e n = + . Aggregating over iV time steps, at t/v = = T, 

Xjv = X(t N ) exp(S N ), (19) 
JV-l 

5jv = e « = fcK V(P,A,A)T + ik K 3 a{p, A, A)W T + o(A;k 4 ), (20) 

n=0 

where ~ N (0, T) and p and a are functions of the parameters determined by 
el°^ and el 1 ^ above, and which are bounded for fixed A > 0. 



9 



For the large wavenumber regime k > h m , for any m > 0, a similar calculation 
to the one in Section 12.21 shows that, under the mean-square stability condition 

X N = o(k p ) Vp>0. 



Following the analysis in Carter fc Giles(2007)| for the deterministic case (p = 
in the present setting), using discrete and continuous Fourier pairs, 



/ir/h 
X n (n) exp(iKhj) d/c, 
-ir/h 



w/h 
-ir/h 

POO 

v(t n ,jh) = t}- X(t n ,K)exp(iKhj) dre, 



2tt 

we can use the decay of X n and X for large k to deduce 

riV „jrr „-,\ 1 



7T//1 

V/ v - v(T,jh) = ± I (X n ( K ) - X(T, k)) exp(iKhj) <1k + o{k) 



j_ 

2tt 



-n/h 
h~ m 

X(T, k) (exp(S'Ar) — 1) exp(inhj) <1k + o(k), 

-h- m 



and together with (|20p we obtain the result upon expanding exp and integrating. 
□ 



3 Convergence tests 

We now test the accuracy and stability of the scheme numerically. 

The computations were conducted with the following set of parameters 
for ((2D, taken from |Bush et al.(2QTT}\ : p = 0.2, /i = 0.081. 



As initial data, we use v(0,x) = S(x — x ) with x = 5. In this case, the 
exact solution to the SPDE can be seen to be 



1 / (x — x — /iT — y/p M> 



" (T ' X,= ^ { l- P )T ^{- '~ ~"2<1- P )t"~ T )■ < 21) 

For the computations, we localise the range of x values to [-16/3,16] 
and set homogeneous Dirichlet boundary conditions. This has been seen to 
introduce negligible localisation error numerically for the above parameters. 

The operators D\ and D2 in ()9]) now have to be interpreted as finite 
difference matrices including the boundary conditions for the first and last 
element. In every timestep, the scheme requires the solution of a tridiagonal 
linear system similar to the one for the heat equation, and therefore has the 
same computational complexity (linear in the number of grid points) as the 
explicit scheme. 



\2 
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The analytical solution u allows us to approximate the mean-square L 2 - 
error using J mesh intervals and N timesteps by 



E(h, kf 



E 



J2(v?(u)-v(Nk,jh;u)) 2 h 



3=0 
M J 



(22) 



m=l j=0 



where the expectation over Brownian paths u> is approximated by the average 
over M samples u m . 

Anticipating applications where the exact solution is unknown, we also 
define error measures based on a fine grid solution / with mesh parameters 
k and h, and a coarse solution c with mesh parameters 4k and 2h, 



e{h, kf 



E 



J/2 

j=0 
M J/2 



(23) 



m=l j=0 



Iterating the refinement, we get a sequence of decreasing grid sizes hi = 
h 2~ l and timesteps k\ = k 4~ l , and denote E\ = E(h t , k t ) the mean-square 
L2-error at level I > and e\ = e(hi, k{). In the following example, ho = 4/3, 
ko = 1/4. The refinement factors are determined by the stability constraint 
of the explicit scheme for k/h 2 and the 0(k, h 2 ) convergence order. 

Note that x does not coincide with a grid point. We apply the Dirac 
initial data to a basis of hat functions to retain second order convergence, 
Pooley et of. (2003)] . 



sec 



Fig. [2] shows the computed values of Ef and e 2 for the explicit scheme. 
6 = a = 0, the drift implicit scheme 9 = 1, a = 0, and the 'Milstein- anti- 
implicit' Crank-Nicolson scheme, 9 = 0.5, a = — 1. The choice of hi and ki is 
within the stability region of all schemes. The results confirm the theoretical 
0(k,h 2 ) convergence, and show that the errors are very similar indeed for all 
schemes. 
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Figure 2: Mean-square error measures E\ and e\ (as explained in the text) 
for the explicit ('expl.', 9 = 0, o = 0), drift implicit ('impl.', 9 = 1, a = 0) 
and Crank-Nicolson-type ('C.N.', # = 0.5, a = —1) Milstein schemes. 

4 An application with locally refined meshes 

4.1 An initial-boundary value problem 

In this section, we consider the initial-boundary value problem on the positive 
half-line, 



i.e., with an absorbing boundary condition at 0. 

We use the same data as in Section [3j and for the numerical tests solve on 
[0, 16] to approximate the positive half-line. The value of 16 was chosen large 
enough that truncation experimentally had negligible impact on the results. 
To illustrate, for T = 5, the standard deviation of each X\ in fl3]) is \fh ~ 2.2, 
so 16 is approximately 5 standard deviations away from their starting point 
Xq = 5. In contrast, the absorbing boundary at is just over 2 standard 
deviations away, which suggests the fraction of absorbed particles (lost mass 
of v) should be in the order of magnitude of 5%. 




(24) 

(25) 
(26) 
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Figure 3: Mean-square error measure e\ (as explained in the text) for the 
explicit ('expl.', 9 = 0, a = 0), drift implicit ('impl.', 9 = 1, a = 0) and 
Crank-Nicolson-type ('C.N.', 9 = 0.5, a = —1) Milstein schemes, for the 
bounded case in comparison with the unbounded case already seen in Fig. [2j 

Fig. El shows that although the estimated error is still asymptotically 
of the same order in this case, it is substantially larger. It is known from 
Krylov fc Lototsky(1998)] that the solution on the half-line is only in H 1 in 
space but does not have L 2 second derivative, however xu xx G L 2 . 

4.2 Local mesh refinement 

To remove the singularity of the computed solution at i = 0, one might 
introduce local coordinate stretching, i.e., a new coordinate y and increasing 
one-to-one function / : [0, oo) — ► [0, oo) with inverse g such that 

v(t,x) = w(tj(x)) & v(t,g(y)) = w(t,y), Vt,x,y>0. 

The SPDE fl2]) in y-coordinates reads 

d^ = (-Ai fog + \ fog} ^ dt + \ (fo gf ^dt-^pfog^ dM t , 

(27) 

and the Milstein finite difference schemes are defined accordingly. 

Conversely, this is closely related to a discretisation of the original SPDE 
on a non-uniform mesh with nodes x n = g(nh). 
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A distinct choice of transformation is y — y/x, because then 

d 2 w d 2 v dv 
= 4 X |- 2 

dy 2 dx 2 dx ' 

and from [Krylov(1994)| the right-hand-side is known to be square-integrable 
in x. This does not imply, however, that w £ H 2 and does not lend itself 
easily to an improved numerical analysis. 

We now investigate the numerical improvement in accuracy for a specific 
application. 



4.3 Application to credit derivatives 



Bush et al. (201 1)) show how the equations (1241) to (|26|) can be used to model 



credit baskets: there, f[2~4"|) describes the evolution of a firm value distribution 
of a large basket of defaultable obligors, where each firm value follows ([3]); 
xq in f[25]) is the firm value at the initial time; fl26|) models default of a firm 
when its value process crosses a default threshold at x — 0. The basket loss, 
i.e. the fraction of firms that have not survived, is then given by 



OG 



1- / u(t,x)dx. (21 



o 

The loss model can be used as the basis for the valuation of basket credit 
derivatives, which are structured by a sequence of regular fee payments by 
the protection buyer, in return for a payment by the protection seller if a 
certain default event occurs. A standardised such product is a collateralized 
debt obligation (CDO) where the payments depend on the losses in a certain 
segment of the basket, measured by attachment points a and detachment 
points d, over a certain time horizon T. The outstanding tranche notional is 
then defined as 

Z t = max(d — L t , 0) — max(a — L t , 0). (29) 

We will consider a maturity T = 5 and a single tranche [a, d] = [0, 0.03]. 

A survey of products and models can be found e.g. in |S chonbucher (2003)] , 
a derivation of pricing formulae in the present model in |Bush et al.(2011)| . 

The main quantities that enter the formulae for tranche spreads are ex- 
pected, discounted (with interest rate r, here 0.042) spread payments, so we 
will be considering here 



p = j2 e- rT nzT^ - z Ti \ 
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Figure 4: Logarithm of variance V\ and mean E[p — Pj_i] of the correction 
from refinement level I — 1 to level Z, for a power grid streching x a . Shown is 
the effect of the change from the original coordinates (equivalent to the case 
a = 1) to the square root stretching (for which a = 1/2). 



where Tj = iq, q = 0.25 (quarterly payments), n = 20. 

We now give results with and without coordinate transformation. For the 
sake of completeness, we provide the approximated loss function in trans- 
formed coordinates 



j-i 

w(t, y)g'(y) dy^l-h^ VjWiVjj: 



where we use the last expression as numerical approximation to the losses. 

Let Pi be an approximation to P with mesh size hi = ho 2~ l and time step 
ki = k 4~ l for I > and ho = 8/5, k = 1/4. We use the 6-a scheme with 
6 = 0.5 and a = — 1 for its unconditional stability. We compute estimators Y/ 
to E[Pj — Pi-i] for I > in order to estimate the contributions of individual 
refinement levels. We do this by averaging p — p_ x over Ni sample paths 
(M t )te[o,T\ (identical for p and P_i, but independent for different Y t ). We 
show E[YJ] and Vi = N{V\Yl] in Fig. HJ where the number of samples, iVj, is 
chosen to make the simulation error negligible. 

The numerical results show that on coarse levels, both the mean and the 
variance of the estimators are much smaller when the computation is done in 
y coordinates with y = y/x instead of the original x coordinates. The smaller 
variance can be exploited by writing, in the spirit of the multi-level Monte 
Carlo method of [Giles(2008)[ |Giles fc Reisinger(2011)] , 



E[P L ] = E 



.1=0 
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for Yi as above for I > and an estimator Y to E[Pq], an d where Vi is the 
variance of Y\ for = 1. 

Using the sum of Y\ as multilevel estimator for E[P], one can optimise 
Ni to give minimal overall computational complexity for a given combined 
variance. As V\ has been reduced on coarser levels due to the grid stretching, 
Ni can be smaller there compared to the uniform grid. Given the complexity 
will be largely determined by the number of samples on the coarsest levels 
(see |Giles fc Reisinger(2011)| ), a decrease in the variance of around 100 on 
those levels (see Fig. H]) allows the number of paths to decrease by a factor 
of 100, which gives significant computational savings. 

5 Conclusions 

We consider implicit variants of the Milstein scheme for a class of SPDEs, 
and show improved stability properties. In particular we find that an 'anti'- 
implicit discretisation of the deterministic part of the Milstein correction 
leads to an unconditionally stable scheme. This is of some importance for stiff 
systems arising from the SPDE discretisation, especially for locally refined 
meshes, where noticable computational savings are observed. 

An important open question is a complete analysis of the numerical ap- 
proximation of initial-boundary value problems for the considered SPDE. It 
might also be interesting to investigate similar ideas in the context of higher 
order expansions of the stochastic integral. 
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